#ifndef MT2_BISECT_H
#define MT2_BISECT_H


/*The user can change the desired precision below, the larger one of the following two definitions is used. Relative precision less than 0.00001 is not guaranteed to be achievable--use with caution*/ 

#define RELATIVE_PRECISION 0.00001 //defined as precision = RELATIVE_PRECISION * scale, where scale = max{Ea, Eb}
#define ABSOLUTE_PRECISION 0.0     //absolute precision for mt2, unused by default


//Reserved for expert
#define MIN_MASS  0.1   //if ma<MINMASS and mb<MINMASS, use massless code
#define ZERO_MASS 0.000 //give massless particles a small mass
#define SCANSTEP 0.1
namespace mt2_bisect
{
class mt2
{  
   public:

      mt2();
      void   mt2_bisect();
      void   mt2_massless();
      void   set_momenta(double *pa0, double *pb0, double* pmiss0);
      void   set_mn(double mn);
      double get_mt2();
      void   print();
      int    nevt;
   private:  

      bool   solved;
      bool   momenta_set;
      double mt2_b;

      int    nsols(double Dsq);
      int    nsols_massless(double Dsq);
      inline int    signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5);
      inline int    signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5);
      int scan_high(double &Deltasq_high);
      int find_high(double &Deltasq_high);
      //data members
      double pax, pay, ma, Ea;
      double pmissx, pmissy;
      double pbx, pby, mb, Eb;
      double mn, mn_unscale;
     
      //auxiliary definitions
      double masq, Easq;
      double mbsq, Ebsq;
      double pmissxsq, pmissysq;
      double mnsq;

      //auxiliary coefficients
      double a1, b1, c1, a2, b2, c2, d1, e1, f1, d2, e2, f2;
      double d11, e11, f12, f10, d21, d20, e21, e20, f22, f21, f20;

      double scale;
      double precision;
};


mt2::mt2()
{
   solved = false;
   momenta_set = false;
   mt2_b  = 0.;
   scale = 1.;
}

double mt2::get_mt2()
{
   if (!momenta_set)
   {
       cout <<" Please set momenta first!" << endl;
       return 0;
   }

   if (!solved) mt2_bisect();
   return mt2_b*scale;
}

void mt2::set_momenta(double* pa0, double* pb0, double* pmiss0)
{
   solved = false;     //reset solved tag when momenta are changed.
   momenta_set = true;

   ma = fabs(pa0[0]);  // mass cannot be negative

   if (ma < ZERO_MASS) ma = ZERO_MASS;

   pax  = pa0[1];
   pay  = pa0[2];
   masq = ma*ma;
   Easq = masq+pax*pax+pay*pay;
   Ea   = sqrt(Easq);

   mb = fabs(pb0[0]);

   if (mb < ZERO_MASS) mb = ZERO_MASS;

   pbx  = pb0[1];
   pby  = pb0[2];
   mbsq = mb*mb;
   Ebsq = mbsq+pbx*pbx+pby*pby;
   Eb   = sqrt(Ebsq);

   pmissx   = pmiss0[1]; pmissy = pmiss0[2];
   pmissxsq = pmissx*pmissx;
   pmissysq = pmissy*pmissy;

// set ma>= mb
   if(masq < mbsq)
   {
      double temp;
      temp = pax;  pax  = pbx;  pbx  = temp;
      temp = pay;  pay  = pby;  pby  = temp;
      temp = Ea;   Ea   = Eb;   Eb   = temp;
      temp = Easq; Easq = Ebsq; Ebsq = temp;
      temp = masq; masq = mbsq; mbsq = temp;
      temp = ma;   ma   = mb;   mb   = temp;
   }
//normalize max{Ea, Eb} to 100
   if (Ea > Eb) scale = Ea/100.;
   else scale = Eb/100.;

   if (sqrt(pmissxsq+pmissysq)/100 > scale) scale = sqrt(pmissxsq+pmissysq)/100;
   //scale = 1;
   double scalesq = scale * scale;
   ma  = ma/scale;
   mb  = mb/scale;
   masq = masq/scalesq;
   mbsq = mbsq/scalesq;
   pax = pax/scale; pay = pay/scale;
   pbx = pbx/scale; pby = pby/scale;
   Ea  = Ea/scale;  Eb = Eb/scale;

   Easq = Easq/scalesq;
   Ebsq = Ebsq/scalesq;
   pmissx = pmissx/scale;
   pmissy = pmissy/scale;
   pmissxsq = pmissxsq/scalesq;
   pmissysq = pmissysq/scalesq;
   mn   = mn_unscale/scale;
   mnsq = mn*mn;

   if (ABSOLUTE_PRECISION > 100.*RELATIVE_PRECISION) precision = ABSOLUTE_PRECISION;
   else precision = 100.*RELATIVE_PRECISION;
}

void mt2::set_mn(double mn0)
{
   solved = false;    //reset solved tag when mn is changed.
   mn_unscale   = fabs(mn0);  //mass cannot be negative
   mn = mn_unscale/scale;
   mnsq = mn*mn;
}

void mt2::print()
{
   cout << " pax = " << pax*scale << ";   pay = " << pay*scale << ";   ma = " << ma*scale <<";"<< endl;
   cout << " pbx = " << pbx*scale << ";   pby = " << pby*scale << ";   mb = " << mb*scale <<";"<< endl;
   cout << " pmissx = " << pmissx*scale << ";   pmissy = " << pmissy*scale <<";"<< endl;
   cout << " mn = " << mn_unscale<<";" << endl;
}

//special case, the visible particle is massless
void mt2::mt2_massless()
{

//rotate so that pay = 0
   double theta,s,c;
   theta = atan(pay/pax);
   s = sin(theta);
   c = cos(theta);

   double pxtemp,pytemp;
   Easq   = pax*pax+pay*pay;
   Ebsq   = pbx*pbx+pby*pby;
   Ea     = sqrt(Easq);
   Eb     = sqrt(Ebsq);

   pxtemp = pax*c+pay*s;
   pax    = pxtemp;
   pay    = 0;
   pxtemp = pbx*c+pby*s;
   pytemp = -s*pbx+c*pby;
   pbx    = pxtemp;
   pby    = pytemp;
   pxtemp = pmissx*c+pmissy*s;
   pytemp = -s*pmissx+c*pmissy;
   pmissx = pxtemp;
   pmissy = pytemp;

   a2  = 1-pbx*pbx/(Ebsq);
   b2  = -pbx*pby/(Ebsq);
   c2  = 1-pby*pby/(Ebsq);

   d21 = (Easq*pbx)/Ebsq;
   d20 = - pmissx +  (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
   e21 = (Easq*pby)/Ebsq;
   e20 = - pmissy +  (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
   f22 = -(Easq*Easq/Ebsq);
   f21 = -2*Easq*(pbx*pmissx + pby*pmissy)/Ebsq;
   f20 = mnsq + pmissxsq + pmissysq - (pbx*pmissx + pby*pmissy)*(pbx*pmissx + pby*pmissy)/Ebsq;

   double Deltasq0    = 0;
   double Deltasq_low, Deltasq_high;
   int    nsols_high, nsols_low;

   Deltasq_low = Deltasq0 + precision;
   nsols_low = nsols_massless(Deltasq_low);

   if(nsols_low > 1)
   {
      mt2_b = (double) sqrt(Deltasq0+mnsq);
      return;
   }

/*
   if( nsols_massless(Deltasq_high) > 0 )
   {
      mt2_b = (double) sqrt(mnsq+Deltasq0);
      return;
   }*/

//look for when both parablos contain origin
   double Deltasq_high1, Deltasq_high2;
   Deltasq_high1 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy;
   Deltasq_high2 = 2*Ea*mn;

   if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high2;
     else Deltasq_high = Deltasq_high1;

   nsols_high=nsols_massless(Deltasq_high);

   int foundhigh;
   if (nsols_high == nsols_low)
   {


      foundhigh=0;

      double minmass, maxmass;
      minmass  = mn ;
      maxmass  = sqrt(mnsq + Deltasq_high);
      for(double mass = minmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
      {
	 Deltasq_high = mass*mass - mnsq;

         nsols_high = nsols_massless(Deltasq_high);
         if(nsols_high>0)
         {
            foundhigh=1;
            Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
            break;
         }
      }
      if(foundhigh==0)
      {

	cout<<"Deltasq_high not found at event " << nevt <<endl;


         mt2_b = (double)sqrt(Deltasq_low+mnsq);
         return;
      }
   }

   if(nsols_high == nsols_low)
   {
      cout << "error: nsols_low=nsols_high=" << nsols_high << endl;
      cout << "Deltasq_high=" << Deltasq_high << endl;
      cout << "Deltasq_low= "<< Deltasq_low << endl;

      mt2_b = sqrt(mnsq + Deltasq_low);
      return;
   }
   double minmass, maxmass;
   minmass = sqrt(Deltasq_low+mnsq);
   maxmass = sqrt(Deltasq_high+mnsq);
   while(maxmass - minmass > precision)
   {
      double Delta_mid, midmass, nsols_mid;
      midmass   = (minmass+maxmass)/2.;
      Delta_mid = midmass * midmass - mnsq;
      nsols_mid = nsols_massless(Delta_mid);
      if(nsols_mid != nsols_low) maxmass = midmass;
      if(nsols_mid == nsols_low) minmass = midmass;
   }
   mt2_b = minmass;
   return;
}

int mt2::nsols_massless(double Dsq)
{
  double delta;
  delta = Dsq/(2*Easq);
  d1    = d11*delta;
  e1    = e11*delta;
  f1    = f12*delta*delta+f10;
  d2    = d21*delta+d20;
  e2    = e21*delta+e20;
  f2    = f22*delta*delta+f21*delta+f20;

  double a,b;
  if (pax > 0) a = Ea/Dsq;
  else         a = -Ea/Dsq;
  if (pax > 0) b = -Dsq/(4*Ea)+mnsq*Ea/Dsq;
  else         b = Dsq/(4*Ea)-mnsq*Ea/Dsq;

  double A4,A3,A2,A1,A0;

  A4 = a*a*a2;
  A3 = 2*a*b2/Ea;
  A2 = (2*a*a2*b+c2+2*a*d2)/(Easq);
  A1 = (2*b*b2+2*e2)/(Easq*Ea);
  A0 = (a2*b*b+2*b*d2+f2)/(Easq*Easq);

  long  double A0sq, A1sq, A2sq, A3sq, A4sq;
  A0sq = A0*A0;
  A1sq = A1*A1;
  A2sq = A2*A2;
  A3sq = A3*A3;
  A4sq = A4*A4;

  long double B3, B2, B1, B0;
  B3 = 4*A4;
  B2 = 3*A3;
  B1 = 2*A2;
  B0 = A1;
  long double C2, C1, C0;
  C2 = -(A2/2 - 3*A3sq/(16*A4));
  C1 = -(3*A1/4. -A2*A3/(8*A4));
  C0 = -A0 + A1*A3/(16*A4);
  long double  D1, D0;
  D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
  D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;

  long double E0;
  E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;

  long  double t1,t2,t3,t4,t5;

//find the coefficients for the leading term in the Sturm sequence
   t1 = A4;
   t2 = A4;
   t3 = C2;
   t4 = D1;
   t5 = E0;

   int nsol;
   nsol = signchange_n(t1,t2,t3,t4,t5)-signchange_p(t1,t2,t3,t4,t5);
   if( nsol < 0 ) nsol=0;

   return nsol;

}

void mt2::mt2_bisect()
{


   solved = true;
   cout.precision(11);

//if masses are very small, use code for massless case.
   if(masq < MIN_MASS && mbsq < MIN_MASS)
   {
      mt2_massless();
      return;
   }


   double Deltasq0;
   Deltasq0 = ma*(ma + 2*mn); //The minimum mass square to have two ellipses

// find the coefficients for the two quadratic equations when Deltasq=Deltasq0.

   a1 = 1-pax*pax/(Easq);
   b1 = -pax*pay/(Easq);
   c1 = 1-pay*pay/(Easq);
   d1 = -pax*(Deltasq0-masq)/(2*Easq);
   e1 = -pay*(Deltasq0-masq)/(2*Easq);
   a2 = 1-pbx*pbx/(Ebsq);
   b2 = -pbx*pby/(Ebsq);
   c2 = 1-pby*pby/(Ebsq);
   d2 = -pmissx+pbx*(Deltasq0-mbsq)/(2*Ebsq)+pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
   e2 = -pmissy+pby*(Deltasq0-mbsq)/(2*Ebsq)+pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
   f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq0-mbsq)/(2*Eb)+
        (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq0-mbsq)/(2*Eb)+
        (pbx*pmissx+pby*pmissy)/Eb)+mnsq;

// find the center of the smaller ellipse
   double x0,y0;
   x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
   y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);


// Does the larger ellipse contain the smaller one?
   double dis=a2*x0*x0+2*b2*x0*y0+c2*y0*y0+2*d2*x0+2*e2*y0+f2;

   if(dis<=0.01)
   {
      mt2_b  = (double) sqrt(mnsq+Deltasq0);
      return;
   }


/* find the coefficients for the two quadratic equations           */
/* coefficients for quadratic terms do not change                  */
/* coefficients for linear and constant terms are polynomials of   */
/*       delta=(Deltasq-m7sq)/(2 E7sq)                             */
   d11 = -pax;
   e11 = -pay;
   f10 = mnsq;
   f12 = -Easq;
   d21 = (Easq*pbx)/Ebsq;
   d20 = ((masq - mbsq)*pbx)/(2.*Ebsq) - pmissx +
         (pbx*(pbx*pmissx + pby*pmissy))/Ebsq;
   e21 = (Easq*pby)/Ebsq;
   e20 = ((masq - mbsq)*pby)/(2.*Ebsq) - pmissy +
         (pby*(pbx*pmissx + pby*pmissy))/Ebsq;
   f22 = -Easq*Easq/Ebsq;
   f21 = (-2*Easq*((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb))/Eb;
   f20 = mnsq + pmissx*pmissx + pmissy*pmissy -
         ((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb)
         *((masq - mbsq)/(2.*Eb) + (pbx*pmissx + pby*pmissy)/Eb);

//Estimate upper bound of mT2
//when Deltasq > Deltasq_high1, the larger encloses the center of the smaller
   double p2x0,p2y0;
   double Deltasq_high1;
   p2x0 = pmissx-x0;
   p2y0 = pmissy-y0;
   Deltasq_high1 = 2*Eb*sqrt(p2x0*p2x0+p2y0*p2y0+mnsq)-2*pbx*p2x0-2*pby*p2y0+mbsq;

//Another estimate, if both ellipses enclose the origin, Deltasq > mT2

   double Deltasq_high2, Deltasq_high21, Deltasq_high22;
   Deltasq_high21 = 2*Eb*sqrt(pmissx*pmissx+pmissy*pmissy+mnsq)-2*pbx*pmissx-2*pby*pmissy+mbsq;
   Deltasq_high22 = 2*Ea*mn+masq;

   if ( Deltasq_high21 < Deltasq_high22 ) Deltasq_high2 = Deltasq_high22;
   else Deltasq_high2 = Deltasq_high21;

//pick the smaller upper bound
   double Deltasq_high;
   if(Deltasq_high1 < Deltasq_high2) Deltasq_high = Deltasq_high1;
   else Deltasq_high = Deltasq_high2;


   double Deltasq_low; //lower bound
   Deltasq_low = Deltasq0;

//number of solutions at Deltasq_low should not be larger than zero
   if( nsols(Deltasq_low) > 0 )
   {
     //cout << "nsolutions(Deltasq_low) > 0"<<endl;
     mt2_b = (double) sqrt(mnsq+Deltasq0);
     return;
   }

   int nsols_high, nsols_low;

   nsols_low  = nsols(Deltasq_low);
   int foundhigh;


//if nsols_high=nsols_low, we missed the region where the two ellipse overlap
//if nsols_high=4, also need a scan because we may find the wrong tangent point.

   nsols_high = nsols(Deltasq_high);

   if(nsols_high == nsols_low || nsols_high == 4)
   {
      //foundhigh = scan_high(Deltasq_high);
      foundhigh = find_high(Deltasq_high);
      if(foundhigh == 0)
      {
 	 cout << "Deltasq_high not found at event " << nevt << endl;
         mt2_b = sqrt( Deltasq_low + mnsq );
         return;
      }

   }

   while(sqrt(Deltasq_high+mnsq) - sqrt(Deltasq_low+mnsq) > precision)
   {
      double Deltasq_mid,nsols_mid;
      //bisect
      Deltasq_mid = (Deltasq_high+Deltasq_low)/2.;
      nsols_mid = nsols(Deltasq_mid);
      // if nsols_mid = 4, rescan for Deltasq_high
      if ( nsols_mid == 4 )
      {
         Deltasq_high = Deltasq_mid;
         //scan_high(Deltasq_high);
         find_high(Deltasq_high);
         continue;
      }


      if(nsols_mid != nsols_low) Deltasq_high = Deltasq_mid;
      if(nsols_mid == nsols_low) Deltasq_low  = Deltasq_mid;
   }
   mt2_b = (double) sqrt( mnsq + Deltasq_high);
   return;
}

int mt2::find_high(double & Deltasq_high)
{
   double x0,y0;
   x0 = (c1*d1-b1*e1)/(b1*b1-a1*c1);
   y0 = (a1*e1-b1*d1)/(b1*b1-a1*c1);
   double Deltasq_low = (mn + ma)*(mn + ma) - mnsq;
   do
   {
      double Deltasq_mid = (Deltasq_high + Deltasq_low)/2.;
      int nsols_mid = nsols(Deltasq_mid);
      if ( nsols_mid == 2 )
      {
         Deltasq_high = Deltasq_mid;
         return 1;
      }
      else if (nsols_mid == 4)
      {
         Deltasq_high = Deltasq_mid;
         continue;
      }
      else if (nsols_mid ==0)
      {
         d1 = -pax*(Deltasq_mid-masq)/(2*Easq);
         e1 = -pay*(Deltasq_mid-masq)/(2*Easq);
         d2 = -pmissx + pbx*(Deltasq_mid - mbsq)/(2*Ebsq)
              + pbx*(pbx*pmissx+pby*pmissy)/(Ebsq);
         e2 = -pmissy + pby*(Deltasq_mid - mbsq)/(2*Ebsq)
              + pby*(pbx*pmissx+pby*pmissy)/(Ebsq);
         f2 = pmissx*pmissx+pmissy*pmissy-((Deltasq_mid-mbsq)/(2*Eb)+
              (pbx*pmissx+pby*pmissy)/Eb)*((Deltasq_mid-mbsq)/(2*Eb)+
              (pbx*pmissx+pby*pmissy)/Eb)+mnsq;
// Does the larger ellipse contain the smaller one?
         double dis = a2*x0*x0 + 2*b2*x0*y0 + c2*y0*y0 + 2*d2*x0 + 2*e2*y0 + f2;
         if (dis < 0) Deltasq_high = Deltasq_mid;
           else Deltasq_low = Deltasq_mid;
      }

   } while ( Deltasq_high - Deltasq_low > 0.001);
   return 0;
}
int mt2::scan_high(double & Deltasq_high)
{
   int foundhigh = 0 ;
   int nsols_high;


   double Deltasq_low;
   double tempmass, maxmass;
   tempmass = mn + ma;
   maxmass  = sqrt(mnsq + Deltasq_high);
   if (nevt == 32334) cout << "Deltasq_high = " << Deltasq_high << endl;
   for(double mass = tempmass + SCANSTEP; mass < maxmass; mass += SCANSTEP)
   {
      Deltasq_high = mass*mass - mnsq;
      nsols_high   = nsols(Deltasq_high);

      if( nsols_high > 0)
      {
	 Deltasq_low = (mass-SCANSTEP)*(mass-SCANSTEP) - mnsq;
         foundhigh   = 1;
         break;
      }
    }
    return foundhigh;
}
int mt2::nsols(  double Dsq)
{
   double delta = (Dsq-masq)/(2*Easq);

//calculate coefficients for the two quadratic equations
   d1 = d11*delta;
   e1 = e11*delta;
   f1 = f12*delta*delta+f10;
   d2 = d21*delta+d20;
   e2 = e21*delta+e20;
   f2 = f22*delta*delta+f21*delta+f20;

//obtain the coefficients for the 4th order equation
//devided by Ea^n to make the variable dimensionless
   long double A4, A3, A2, A1, A0;

   A4 =
   -4*a2*b1*b2*c1 + 4*a1*b2*b2*c1 +a2*a2*c1*c1 +
   4*a2*b1*b1*c2 - 4*a1*b1*b2*c2 - 2*a1*a2*c1*c2 +
   a1*a1*c2*c2;

   A3 =
     (-4*a2*b2*c1*d1 + 8*a2*b1*c2*d1 - 4*a1*b2*c2*d1 - 4*a2*b1*c1*d2 +
   8*a1*b2*c1*d2 - 4*a1*b1*c2*d2 - 8*a2*b1*b2*e1 + 8*a1*b2*b2*e1 +
   4*a2*a2*c1*e1 - 4*a1*a2*c2*e1 + 8*a2*b1*b1*e2 - 8*a1*b1*b2*e2 -
     4*a1*a2*c1*e2 + 4*a1*a1*c2*e2)/Ea;


   A2 =
     (4*a2*c2*d1*d1 - 4*a2*c1*d1*d2 - 4*a1*c2*d1*d2 + 4*a1*c1*d2*d2 -
   8*a2*b2*d1*e1 - 8*a2*b1*d2*e1 + 16*a1*b2*d2*e1 +
   4*a2*a2*e1*e1 + 16*a2*b1*d1*e2 - 8*a1*b2*d1*e2 -
   8*a1*b1*d2*e2 - 8*a1*a2*e1*e2 + 4*a1*a1*e2*e2 - 4*a2*b1*b2*f1 +
   4*a1*b2*b2*f1 + 2*a2*a2*c1*f1 - 2*a1*a2*c2*f1 +
     4*a2*b1*b1*f2 - 4*a1*b1*b2*f2 - 2*a1*a2*c1*f2 + 2*a1*a1*c2*f2)/Easq;

   A1 =
     (-8*a2*d1*d2*e1 + 8*a1*d2*d2*e1 + 8*a2*d1*d1*e2 - 8*a1*d1*d2*e2 -
   4*a2*b2*d1*f1 - 4*a2*b1*d2*f1 + 8*a1*b2*d2*f1 + 4*a2*a2*e1*f1 -
   4*a1*a2*e2*f1 + 8*a2*b1*d1*f2 - 4*a1*b2*d1*f2 - 4*a1*b1*d2*f2 -
     4*a1*a2*e1*f2 + 4*a1*a1*e2*f2)/(Easq*Ea);

   A0 =
     (-4*a2*d1*d2*f1 + 4*a1*d2*d2*f1 + a2*a2*f1*f1 +
   4*a2*d1*d1*f2 - 4*a1*d1*d2*f2 - 2*a1*a2*f1*f2 +
     a1*a1*f2*f2)/(Easq*Easq);

   long  double A0sq, A1sq, A2sq, A3sq, A4sq;
   A0sq = A0*A0;
   A1sq = A1*A1;
   A2sq = A2*A2;
   A3sq = A3*A3;
   A4sq = A4*A4;

   long double B3, B2, B1, B0;
   B3 = 4*A4;
   B2 = 3*A3;
   B1 = 2*A2;
   B0 = A1;

   long double C2, C1, C0;
   C2 = -(A2/2 - 3*A3sq/(16*A4));
   C1 = -(3*A1/4. -A2*A3/(8*A4));
   C0 = -A0 + A1*A3/(16*A4);

   long double D1, D0;
   D1 = -B1 - (B3*C1*C1/C2 - B3*C0 -B2*C1)/C2;
   D0 = -B0 - B3 *C0 *C1/(C2*C2)+ B2*C0/C2;

   long double E0;
   E0 = -C0 - C2*D0*D0/(D1*D1) + C1*D0/D1;

   long  double t1,t2,t3,t4,t5;
//find the coefficients for the leading term in the Sturm sequence
   t1 = A4;
   t2 = A4;
   t3 = C2;
   t4 = D1;
   t5 = E0;


//The number of solutions depends on diffence of number of sign changes for x->Inf and x->-Inf
   int nsol;
   nsol = signchange_n(t1,t2,t3,t4,t5) - signchange_p(t1,t2,t3,t4,t5);

//Cannot have negative number of solutions, must be roundoff effect
   if (nsol < 0) nsol = 0;

   return nsol;

}

inline int mt2::signchange_n( long double t1, long double t2, long double t3, long double t4, long double t5)
{
   int nsc;
   nsc=0;
   if(t1*t2>0) nsc++;
   if(t2*t3>0) nsc++;
   if(t3*t4>0) nsc++;
   if(t4*t5>0) nsc++;
   return nsc;
}
inline int mt2::signchange_p( long double t1, long double t2, long double t3, long double t4, long double t5)
{
   int nsc;
   nsc=0;
   if(t1*t2<0) nsc++;
   if(t2*t3<0) nsc++;
   if(t3*t4<0) nsc++;
   if(t4*t5<0) nsc++;
   return nsc;
}



}//end namespace mt2_bisect

#endif
